Is BCS theory exact in the thermodynamic limit? 
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The exact ground state of the reduced BCS Hamiltonian with a non-degenerate single-particle 
spectrum is investigated numerically for large system sizes (up to 4000 levels) and compared with 
the BCS ansatz. While generic quantities such as the total energy or the superconducting order 
parameter are found to converge to the BCS values already for relatively small system sizes, other 
quantities, in particular the ground state fidelity susceptibility and a pair correlation function in- 
volving HOMO and LUMO levels, do not agree with BCS theory, even for large system sizes where 
they have visibly converged to well-defined limits. The discrepancies are particularly large in the 
intermediate coupling range. We conclude that the BCS wave function does not represent the exact 
ground state of the reduced BCS Hamiltonian. 



The highly successful microscopic theory of supercon- 
ductivity by Bardeen, Cooper and Schrieffer [1] involves 
two fundamental elements, on the one hand the so-called 
reduced BCS Hamiltonian, where only scattering pro- 
cesses between zero-momentum pairs of electrons are 
taken into account, on the other hand a variational ansatz 
for the ground state of this Hamiltonian, a coherent su- 
perposition of products of pair wave functions. In this 
paper we address the question whether the BCS ansatz 
is the exact ground state of the reduced BCS Hamilto- 
nian in the limit of an infinitely large system size. This 
problem has already been raised in the original BCS pa- 
per, where the cautious answer "may be even correct 
in the limit of a large number of particles" was given. 
Shortly after, Anderson argued that BCS theory should 
be "nearly valid" because in the limit of large numbers 
"quantum fluctuations die out" [2]. Later explicit cal- 
culations showed that indeed the ground state energy, 
level occupation and the free energy were exactly pre- 
dicted by BCS theory in the thermodynamic limit [3-6]. 
Moreover, for a specific single-particle spectrum ("step 
model"), Mattis and Lieb concluded that the BCS wave 
function was exact in this limit [7]. One is therefore 
tempted to answer with an unconditional yes to the ques- 
tion of asymptotic validity of BCS theory. Below we 
will add further support to this conclusion by showing 
(numerically) that the superconducting order parameter 
tends to the BCS limit for large system sizes. However, 
we also find subtle correlations, such as the ground state 
fidelity susceptibility and a "HOMO-LUMO" pair corre- 
lation function, for which the BCS predictions differ from 
the (numerically) exact result, even for very large system 
sizes. 
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where c\, a and c va are, respectively, fermionic creation 
and annihilation operators for electrons in level v with 
spin a. This Hamiltonian has electron-hole symmetry, 
and therefore the chemical potential vanishes if the par- 
ticle number N equals L (half filling) , the case considered 
in this paper. 

The conventional BCS ground state is defined as 
I^bcs) = ILK + ^4; c If)l°)) where u " = [(&» + 
e v )/{2E u )] 1 l 2 , v v = [{E v - t v )/{2E u )] l l 2 with E„ = 
(e 2 + A 2 ) 1 / 2 . The gap parameter A is determined by 
minimizing the energy expectation value. For the present 
model we obtain A = [2 sinh(l/g)] _1 in the limit L^oo. 
The BCS state can also be written as 
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The exact eigenstates of the Hamiltonian (1) were first 
derived by Richardson [8] , who used it as a model for de- 
scribing nucleons coupled by pairing forces. It was proven 
later [10] that this Hamiltonian is part of a larger family 
of integrable models for which eigenstates and eigenval- 
ues were found by Gaudin [11] in the framework of the 
Bethe ansatz. The eigenstates for N — 2M particles have 
the fornr 
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Our analysis is based on the exact treatment of the 
reduced BCS Hamiltonian by Richardson [8, 9], who 
also proposed a convenient single-particle spectrum, e„ = 
-W/2 +™{v- 1/2) with v = 1, . . . , L. This model of 
equally spaced levels has no level degeneracy and may 
be associated with a system of chiral fermions running 
around a ring in only one direction. We will use the 
width W as unit of energy, W — 1. The reduced BCS 



where the "rapidities" Ai satisfy the Richardson equations 
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These equations can be solved numerically by scanning 
from g — up to some finite value g > 0. Analyti- 
cal insight has been provided by Gaudin [12, 13] in the 



continuum limit (L — >• oo), using an analogy to electro- 
statics. In this limit the rapidities for the ground state 
are found to lie on a well defined curve in the complex 
plane. For equally spaced energy levels this result was 
used to show [13] that the BCS equations for the gap, 
the chemical potential and the ground state energy are 
reproduced in the thermodynamic limit. The low energy 
excitations have also been obtained by solving Eq. (4) 
analytically in the strong coupling limit g — > oo, where 
all (ground state) rapidities are complex and diverge [14]. 
The numerical calculations are greatly simplified by 



introducing the variables A(e„) 



E, 



satisfy the "substituted Bethe equations" [15] 
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This allows us to study much larger system sizes than 
by solving the Richardson equations (4). However some 
particular quantities can only be written in terms of the 
rapidities. In this case, once Eqs. (5) are solved we follow 
the procedure outlined in [16] to extract the rapidities. 

The ground state energy E (g) can be directly calcu- 
lated from the solution of the substituted Bethe equa- 
tions, because it depends on the variables A(e„) through 
the simple relation Eq = ^2 V 2e u h{t v ) — Lg/A. We have 
calculated Eq for various system sizes L and observed 
that it converges rapidly towards the thermodynamic 
limit of the BCS result, E^ cs = -(£/4)coth(l/p). 

A more interesting quantity is the ground state fidelity 
susceptibility, which has been introduced to characterize 
quantum phase transitions [17, 18], and which can also 
be used to describe crossover phenomena [19, 20]. The 
overlap F(g, Sg) = (^> g \^> g+ s g ) between the ground states 
obtained for two different model parameters, in our case 
two different coupling constants, is the ground state fi- 
delity. It is convenient to consider the case of infinitely 
close parameters and to define the fidelity susceptibility 
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In BCS theory the fidelity susceptibility can be obtained 
analytically. For the case studied here we obtain (for 
L -> oo) 
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In the Bethe ansatz framework the fidelity F(g,Sg) = 
(fygl'fyg+Sg} is given by the determinant of an L x L ma- 
trix [21], from which the fidelity susceptibility is calcu- 
lated using Eq. (6). Fig. 1 shows the exact results 
obtained in this way for different system sizes in com- 
parison with the BCS result for L — >• oo. We observe a 
rapid convergence to a limiting curve for L rs 500. For 



larger systems finite size effects are negligible (except for 
very small couplings) and we can thus assume to have 
reached the thermodynamic limit. Clearly the limiting 
curve disagrees with the BCS result. The difference is 
maximal around g « 0.26 which is also the value where 
both curves have a maximum. 




FIG. 1. Ground state fidelity susceptibility as a function 
of coupling strength g. Symbols represent solutions of the 
substituted Bethe equations (5) for various values of L, the 
dashed line stands for the BCS result for L — > oo. 

The ground state fidelity susceptibility can be repre- 
sented with respect to the eigenstates {\^ n (g))} of the 
Hamiltonian with coupling constant g, by using ordinary 
perturbation theory in powers of Sg. One finds 



x(g) 






|(tto(g»M|ttn(g))| S 

[EQ{g)-E n (g)Y 



(8) 



The energy eigenvalues E n (g) are expected to converge 
to the BCS values for L — > oo and therefore the dis- 
crepancy between exact and BCS results must be related 
to pair-pair correlations. Moreover the effect should be 
largest for small excitation energies because of the de- 
nominator in Eq. (8). Therefore we have also evaluated 
the correlation function 
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choosing /i to be the highest occupied level (HOMO) and 
v the lowest unoccupied level (LUMO), in the case g = 0. 
This function is easily calculated for the BCS ground 
state, where it is given by S(v, (j) = (LA) 2 /[1 + (2LA) 2 }. 
A represents the gap parameter for L levels (and N = L 
electrons), which vanishes for g < g c (L). For large L 
an excellent approximation for the (BCS) critical cou- 
pling strength is g c (L) « [log (2L) + 7] -1 , where 7 is 
Euler's constant (7 = 0.5772157...). For L ->• 00 the 
BCS expression tends to 1/4 for any finite coupling g, 
while it vanishes for g = 0. In the exact solution shown 
in Fig. 2 g c (L) coincides approximately with the region 
where S increases most rapidly. For smaller values of g 
it does not make sense to compare BCS (S = 0) with the 
exact behavior, which can be calculated analytically for 



g — > (and L finite), namely S ~ g/2, in agreement with 
the asymptotic behavior of Fig. 2. For large couplings 
the exact results converge rapidly to the BCS prediction, 
however the exact correlation function exhibits a peak 
which does not decrease with increasing system size, in 
sharp contrast to the BCS prediction. Once again the 
largest discrepancy appears within the crossover region 
between weak and strong coupling. An educated guess 
for the limit L — > oo is a step at g = 0, with a larger 
size than BCS, followed by a smooth decrease towards 
the asymptotic strong-coupling limit 1/4. 




FIG. 2. HOMO-LUMO pair correlation function. Symbols 
represent the exact solution, while the dashed lines stand for 
the BCS result for various L. 

For general labels p, v Eq. (9) is essentially the reduced 
density matrix pi(y,p) := ("Jol^^l^o) introduced by 
Yang [22]. For the BCS ground state we get (for L — > oo) 
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According to Yang there is off-diagonal long-range or- 
der if (and only if) the reduced density matrix has an 
eigenvalue of the order of the particle number, i.e., of or- 
der L in the present case. For the BCS ground state we 
have found that the eigenvalues Ui of p2 are solutions of 
the equation 
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Most solutions for ui are such that the second term 
changes sign when summing from v = 1 to v = L, but 
when this does not happen the solution has to be of order 
L, lo = ah. In the present case we find 
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We have also analyzed the reduced density matrix for 
the exact ground state. The matrix elements can be cal- 
culated as sums of certain determinants which depend 
on the rapidities [23]. Fig. 3 shows that the largest 



eigenvalue is a linear function of L already for rather 
small system sizes and for small coupling strengths. The 
proportionality constant, which is directly related to the 
superconducting order parameter, agrees perfectly well 
with the BCS result. A similar conclusion was recently 
reached by a somewhat different approach [24]. 




FIG. 3. Order parameter a — \im.L-Hx(w ma ,x/L), where u>max 
is the largest eigenvalue of the reduced density matrix p2- The 
exact values (dots) were extracted from the L-dependence of 
Wmax for a few selected coupling constants (inset). The full 
line represents the BCS result. 

In the pseudospin language introduced by Ander- 
son [2] the pair correlation function (9) is nothing but 
(^ o\S£ S~\*f! q) . It is therefore natural to study also the 
quantity (\I/o|S„ • S^l^o)) which can be obtained from the 
solutions of a set of modified substituted Bethe equa- 
tions. In fact, for p ^ v we can write this correlation 
function as the derivative of the expectation value of a 
certain R operator with respect to e u . The R opera- 
tors commute among themselves and with the Hamilto- 
nian (1) and therefore one can use the Hellman-Feynman 
theorem for replacing the derivative of the expectation 
value by the expectation value of the derivative. But the 
derivative of R operators yields essentially scalar prod- 
ucts S„ • S M . Surprisingly, as shown in Fig. 4, this cor- 
relation function appears to tend to the BCS result (for 
g > g c (L)), in contrast to the pair-pair correlation func- 
tion of Fig. 2. We speculate that this difference is related 
to the invariance of S„ • S^ with respect to a rotation in 
pseudospin space. 

We have seen that in the thermodynamic limit many 
exact quantities converge to the values predicted by BCS 
theory. One may wonder whether the remaining discrep- 
ancies, namely those found for both the fidelity suscepti- 
bility and the HOMO-LUMO pair correlations, also dis- 
appear if, instead of the conventional BCS ansatz, we use 



(12) the number-projected state 
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where 



the operator B< has been defined in Eq. (2). This state 
indeed resembles the Richardson solution (3), but, as 
emphasized by Combescot and collaborators [25], in the 
BCS state all pair operators are equal (B'), while they 
are all different in the Richardson solution (Bj). 
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FIG. 4. HOMO-LUMO pseudospin correlation function. 
Symbols represent solutions of substituted Bethe equations 
for derivatives of A(e„), dashed lines depict the BCS predic- 
tions, which tend to the limiting constant 1/4 for L — s> oo. 



It is not an easy task to deal with the number-projected 
BCS ansatz for large systems and therefore we have lim- 
ited ourselves to modest sizes. The results for the fidelity 
susceptibility, shown in Fig. 5 for both conventional and 
projected BCS states, indicate a clear convergence for 
not too small couplings, in contrast to the results dis- 
played in Fig. 1. Therefore the discrepancies cannot be 
removed by replacing the conventional BCS ansatz by 
the number-projected state. 
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FIG. 5. Fidelity susceptibility according to BCS theory. Sym- 
bols represent the results obtained for the number-projected 
BCS ground state for two values of L, while the thermody- 
namic limit for the conventional BCS ansatz is given by the 
dashed line. 

In conclusion, we have studied the reduced BCS Hamil- 
tonian using both the exact ground state and the BCS 
ansatz. We have confirmed that in the thermodynamic 
limit some generic quantities, calculated for the exact 
ground state, converge rapidly towards the values pre- 
dicted by BCS theory with increasing system size. How- 
ever there exist subtle quantities, such as the fidelity sus- 
ceptibility and the HOMO-LUMO pair correlation func- 
tion, where discrepancies remain. The differences are es- 
pecially pronounced in the crossover region from weak to 
strong coupling. In this sense BCS theory is not exact in 



all respects in the thermodynamic limit. Unfortunately, 
it is far from obvious how to measure the fidelity suscep- 
tibility or the HOMO-LUMO pair correlation function 
for an actual material. A way out of this dilemma would 
be to proceed to experimentally accessible dynamic cor- 
relation functions, where the discrepancies between ex- 
act and mean-field predictions may be even more pro- 
nounced. 
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